{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "498d70fa",
   "metadata": {},
   "source": [
    "# 线性回归（二）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "210f3578",
   "metadata": {},
   "source": [
    "<font color=blue size=4>1.课堂实验任务</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "226badd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70d2906a",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">1) 使用pandas库的read_csv()函数(可以参考[pandas的官方文档](https://pandas.pydata.org/docs/reference/api/pandas.read_csv.html))将训练数据集'train.csv'和测试数据集'test.csv'载入到Dataframe对象中。</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f59640f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here\n",
    "\n",
    "#读取数据集\n",
    "train_frame = pd.read_csv('train.csv')\n",
    "test_frame = pd.read_csv('test.csv')\n",
    "\n",
    "#转化成numpy矩阵\n",
    "train = np.array(train_frame)\n",
    "test = np.array(test_frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bccb84d",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">2) 假设模型为一元线性回归模型$\\hat{y}=wx+b$, 损失函数为$l(w,b)=\\frac{1}{2}\\sum_{i=1}^m(\\hat{y}^{(i)}-y^{(i)})^2$, 其中$\\hat{y}^{(i)}$表示第$i$个样本的预测值，$y^{(i)}$表示第$i$个样本的实际标签值, $m$为训练集中样本的个数。求出使得损失函数最小化的参数$w$和$b$。</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e011976d",
   "metadata": {},
   "source": [
    "方法① \n",
    "\n",
    "将$l(w,b)$分别对$w$和$b$求导，得到\n",
    "$$\n",
    "\\frac{\\partial l(w,b)}{\\partial w}=w\\sum_{i=1}^m x_i^2 -\\sum_{i=1}^m (y_i-b)x_i,\n",
    "$$\n",
    "$$\n",
    "\\frac{\\partial l(w,b)}{\\partial b}=mb -\\sum_{i=1}^m (y_i-wx_i),\n",
    "$$\n",
    "令上述两式为零即可得到$w$和$b$的解析解：\n",
    "$$\n",
    "w=\\frac{\\sum_{i=1}^m y_i (x_i-\\bar{x})}{\\sum_{i=1}^m x_i^2-\\frac{1}{m}(\\sum_{i=1}^m x_i)^2},\n",
    "$$\n",
    "$$\n",
    "b=\\frac{1}{m}\\sum_{i=1}^m(y_i-wx_i),\n",
    "$$\n",
    "其中$\\bar{x}=\\frac{1}{m}\\sum_{i=1}^m x_i$为$x$的均值。\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a04bc78",
   "metadata": {},
   "source": [
    "方法② 梯度下降法。手动实现梯度下降法(不使用机器学习框架，如PyTorch、TensorFlow等)来进行模型的训练。算法步骤如下：1.初始化模型参数$w$和$b$的值；2.在负梯度的方向上更新参数(批量梯度下降、小批量随机梯度下降或者随机梯度下降均可)，并不断迭代这一步骤，更新公式(以小批量随机梯度下降为例)可以写成：$$w\\gets w-\\frac{\\eta}{\\left|B\\right|}\\sum_{i\\in{B}}x^{(i)}(wx^{(i)}+b-y^{(i)})$$, 和$$b\\gets b-\\frac{\\eta}{\\left|B\\right|}\\sum_{i\\in{B}}(wx^{(i)}+b-y^{(i)})$$， 其中$\\eta$表示学习率,$B$表示每次迭代中随机抽样的小批量，$\\left|B\\right|$则表示$B$中的样本数量。3. 终止条件为迭代次数达到某一上限或者参数更新的幅度小于某个阈值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bad649ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "func_1: omega = 3.041479, b = 4.906074\n",
      "func_2: omega = 3.088664, b = 4.833292\n"
     ]
    }
   ],
   "source": [
    "# 方法①\n",
    "def func_1():\n",
    "    x_average = np.mean(train[:, 0])\n",
    "    omega = np.sum((train[:, 0] - x_average) * train[:, 1]) / (np.sum(train[:, 0] ** 2) - np.sum(train[:, 0]) ** 2 / len(train))\n",
    "    b = np.mean(train[:, 1] - omega * train[:, 0])\n",
    "    return omega, b\n",
    "\n",
    "# 方法②\n",
    "# 将数据集分成k个mini_batch\n",
    "def func_2(omega_init, b_init, k, learning_rate, threshold):\n",
    "    def get_mini_batches(train, k):\n",
    "        np.random.shuffle(train)\n",
    "        mini_batches = [train[i:i + k] for i in range(0, len(train), k)]\n",
    "        return mini_batches\n",
    "    \n",
    "    mini_batches = get_mini_batches(train, k)\n",
    "\n",
    "    def gradient_descent(omega, b, mini_batches):\n",
    "        # omega的更新公式\n",
    "        def omega_gd(omega, b, mini_batch):\n",
    "            total = 0\n",
    "            for x, y in mini_batch:\n",
    "                total += x * (omega * x + b - y)\n",
    "            return omega - learning_rate * total / len(mini_batch)\n",
    "\n",
    "        # b的更新公式\n",
    "        def b_gd(omega, b, mini_batch):\n",
    "            total = 0\n",
    "            for x, y in mini_batch:\n",
    "                total += omega * x + b - y\n",
    "            return b - learning_rate * total / len(mini_batch)\n",
    "        \n",
    "        while True:\n",
    "            random_index = np.random.randint(0, len(mini_batches))\n",
    "            mini_batch = mini_batches[random_index]\n",
    "            omega_new, b_new = omega_gd(omega, b, mini_batch), b_gd(omega, b, mini_batch)\n",
    "            # 终止条件为参数更新的幅度小于阈值threshold\n",
    "            if abs(omega_new - omega) < threshold and abs(b_new - b) < threshold:\n",
    "                break\n",
    "            omega, b = omega_new, b_new\n",
    "        return omega, b\n",
    "    return gradient_descent(omega_init, b_init, mini_batches)\n",
    "\n",
    "omega_1, b_1 = func_1()\n",
    "print('func_1: omega = %f, b = %f' % (omega_1, b_1))\n",
    "\n",
    "#omega和b为任意初始值，k取10，学习率取0.01，阈值取0.0001\n",
    "k = 10\n",
    "learning_rate = 0.01\n",
    "threshold = 0.0001\n",
    "\n",
    "omega_2, b_2 = func_2(1, 1, k, learning_rate, threshold)\n",
    "print('func_2: omega = %f, b = %f' % (omega_2, b_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "844f2a87",
   "metadata": {},
   "source": [
    "方法③ \n",
    "\n",
    "用矩阵表示，假设数据集有$m$个样本，特征有$n$维$。X=\\left[ \\begin{matrix} x_{11} & x_{12} & \\cdots & x_{1n} & 1 \\\\\n",
    "                         x_{21} & x_{22} & \\cdots & x_{2n} & 1 \\\\\n",
    "                         \\vdots & \\vdots &      & \\vdots & \\vdots \\\\\n",
    "                         x_{m1} & x_{m2} & \\cdots & x_{mn} & 1 \\end{matrix} \\right]$,\n",
    "        实际标签$Y=\\left[ \\begin{matrix} y_{1} \\\\\n",
    "                         y_{2} \\\\\n",
    "                         \\vdots \\\\\n",
    "                         y_{m}\\end{matrix} \\right]$,\n",
    "        参数$B=\\left[ \\begin{matrix} w_{1} \\\\\n",
    "                         w_{2} \\\\\n",
    "                         \\vdots \\\\\n",
    "                         w_{n} \\\\\n",
    "                         b\\end{matrix} \\right]$，则解析解为$B^*=(X^T X)^{-1}X^T Y$。推导过程可参考[这篇文章](https://zhuanlan.zhihu.com/p/74157986)。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb16c177",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">3) 使用求解出来的线性回归模型对测试数据集'test.csv'进行预测，输出可视化结果（比如用seaborn或者matplotlib等可视化库来画出测试数据的散点图以及训练好的模型函数图像）。</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b19a5413",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0L0lEQVR4nO3de3iU9Z3//9dnBojhkAEShKQkJlrAqoGlghWx1EOVol/WbtYWdat8235bbUXrYa0Ea2vdnwR6WVfW43pYi63W7Bbw5/fXWsSCIouyQUBAWxEbCEooRnESkhgg8/n9kWTCZGaSzOSee07Px3Xlupz7vifz8fa+nFfen5Ox1loBAAC4xJPsBgAAgOxC+AAAAK4ifAAAAFcRPgAAgKsIHwAAwFWEDwAA4CrCBwAAcBXhAwAAuGpQshvQUyAQ0P79+zVixAgZY5LdHAAA0A/WWjU1NamoqEgeT++1jZQLH/v371dxcXGymwEAAOKwb98+jR8/vtdrUi58jBgxQlJH4/Py8pLcGgAA0B+NjY0qLi4Ofo/3JuXCR1dXS15eHuEDAIA0058hEww4BQAAriJ8AAAAVxE+AACAqwgfAADAVYQPAADgKsIHAABwFeEDAAC4ivABAABcRfgAAACuInwAAABXET4AAEhx9f5WbXy/QfX+1mQ3xREpt7cLAADoVl1Tp8qVOxSwksdIVRXlmje9JNnNGhAqHwAApKh6f2sweEhSwEqLVu5M+wpIzOFj/fr1mjt3roqKimSM0fPPPx/12muvvVbGGN1///0DaCIAANmptqE5GDy6tFurPQ0tyWmQQ2IOH83NzZoyZYoefPDBXq97/vnntWnTJhUVFcXdOAAAsllZwTB5euxQ7zVGpQVDk9Mgh8Q85mPOnDmaM2dOr9d8+OGHWrBggVavXq1LL7007sYBAJDNCn25qqoo16KVO9VurbzGaHHFGSr05Sa7aQPi+IDTQCCgq6++WrfddptOP/30Pq9va2tTW1tb8HVjY6PTTQIAIG3Nm16iWRPHaE9Di0oLhg4seBw6JC1YIF16qXTVVc41MkaODzhdunSpBg0apBtvvLFf11dVVcnn8wV/iouLnW4SAABprdCXqxmn5A8seDz3nDR6tOr/72ptvHtZUgetOlr5ePPNN7Vs2TJt2bJFxpi+3yCpsrJSt9xyS/B1Y2MjAQQAAKcEAlJ5ufTOO6qefJEqZ9+ggMcjz5K1SZu262jl47XXXtPBgwdVUlKiQYMGadCgQdq7d69uvfVWlZaWRnxPTk6O8vLyQn4AAIAD/vIXyeuV3nlH9SPyg8FDSu60XUcrH1dffbW++tWvhhybPXu2rr76an3729928qMAAEBvFi2SqqqCL/971LRg8OjSNW3X7QGsMYePw4cPa/fu3cHXtbW12rZtm0aPHq2SkhLl5+eHXD948GCNGzdOkyZNGnhrAQBA71papGHDQg5do+V69tA39LnAWpnj8keypu3G3O2yefNmTZ06VVOnTpUk3XLLLZo6dap++tOfOt44AAAQg9Wrw4JHgT7Sr3WN2pty9Q/jy+XtHJOZzGm7xlpr+77MPY2NjfL5fPL7/Yz/AACgP6yVLr5Yevnl4KFf61u6Rr8Ovv7kE2nUqI4l2x2ZtttDLN/fbCwHAEA627dPKgmdsTJDG/WGZkiSbrxRWras+1yhLzfpi5QRPgAASFdlZdKePcGXrTpBPvl1VEMkSTt3Sv1Y79N17GoLAEC6aWmRjAkJHjfrPg1Vq45qiMrLpfb21AweEuEDAID0smxZ2KDSmdqg+3WzJGnFCmn7dsmTwt/wdLsAAJAuIqwebhSQ1HG8sVEaMcLlNsUhhXMRAADZqd7fqo3vN3SvPrpxY1jwWKYbZWQlGRUXd0x4SYfgIVH5AAAgpVTX1Kly5Q4FrOQxUtUflmne9jUh14zRQTVojCTpD3+Q5sxJRkvjR/gAACBF1Ptbg8FD6tx/ZfYCzardosKmjyWps9rReT4QsScm5dHtAgBAiqhtaA4Gjy7tHq/2jCzSHfp/gsHj+us7ulnSMXhIVD4AAEgZZQXD5AkEQjaA8wbadd6nG3VEwyV1rCk2fnyyWugMKh8AAKSCt95S4cihqlr9gLyBdkkdwePg6r/TkaaO4LH/01bVtR03EDVNUfkAACDZjus/mbd9jWbVbtHlI5/Rpk8vUHtTrp5+Whpyap1mLjluIGpFueZNL+nll6YuwgcAAMlibcTVwIqaGqSmjn8OBKQDja3B4CF1DkRduVOzJo5J+j4t8aDbBQCAZPj+98OCx591anBQadfaHcZEGYhqrfY0tLjVWkdR+QAAwG0RpqmM1sc6pNGSpA8/lIqKus+VFQyTxygkgHiNUWnB0ES3NCGofAAA4Ja9e6MskW6DwcPa0OAhSYW+XFVVlMvb+V6vMVpccUZadrlIVD4AAHBHhNBxjxbpJ7pHkvTww9IPfhD97fOml2jWxDHa09Ci0oKhaRs8JMIHAACJF6XaIUneEa1av7VZZQXDJPUeKAp9uWkdOrrQ7QIAQKI89FCvwWP45DqVXL9WVz2+STOXrFV1TZ3bLUwKwgcAAIlgjLRgQcihyXorGDy2727VmEvCp8+m+wJi/UH4AADASYcPR6127NBkSR2DSpuUWdNnY0H4AADAKYWF0ogRIYfW6KvBasfDD3cED6l7+uzx0nn6bCwYcAoASAn1/lbVNnQMvEzLQZURqh2DdFTtnV+1tkeVo2v67KKVO9VubdpPn40F4QMAkHTVNXWqXJmm+5asWiVVVIQd7qp2DB4sHTkS+a2ZNH02FnS7AACSqt7fGgweUpoNvDQmLHjM16+CwaOxMXrw6FLoy9WMU/KzJnhIVD4AAEnW274lKfuFHAhIXm/Y4a7QIYV3s6AblQ8AQFKl3cDLs8/uNXg8/jjBoy9UPgAASZVWAy8jDCodo4Nq0BhJhI7+InwAAJIu5QdevvmmNG1a2GG6WeJDtwsAICWk7MBLY8KCxwNaEAwen31G8IgVlQ8AAKLpZV8WidARLyofAAD0dM01vQaPF14geAwElQ8AAI4XIXR8SW/of/QlSYQOJxA+AACQpI8+kk48Meww3SzOI3wAABCh2nFYwzRChyVJx45FXNoDcSJ8AABSkmsbzUUIHl4dU0AdaYNqh/MYcAoASDnVNXWauWStrnp8k2YuWavqmjrnP+TRR6MOKg3Iqz/9ieCRKFQ+AAApJdpGc7MmjnGuAhIhdHxXT+g/9F1JhI5EI3wAAFJKQjeaO3pUGjIk7HDXoNIJE6Rdu+L/9a51FaU5wgcAIKV0bTR3fABxZKO5wYM7Ro720BU8AoGIBZF+q66pC1ZsPEaqqijXvOkl8f/CDMaYDwBASunaaM7bmQQc2WjOmLDgka+GYPCwdmDBI1pXUb2/Nf5fmsGofAAAUo5jG839939L554bdrgrdGzbJk2ZMoCGdkpoV1EGInwAAFJSoS934NWOHp7Qd/U9PSHJ2UGlCesqylB0uwAAMk+UKbTf0xO69FLnZ7MkpKsog1H5AABkjssu69j1rYfjx3YkimNdRVmA8AEAyAwRqh1TtE3b1TGow421OwbcVZQl6HYBAKS3Dz+M2s2yXVO0Z48zwaPe36qN7zcwg8UBVD4AAOkrQuioU7FOUsdy7E5VO1jDw1lUPgAA6SlC8PCoXSepTjfc4FzwSPU1PNKxIkPlAwCQXu69V7rttrDDiRpUmspreKRrRYbKBwAgYRz/q9yYsOBxhX6b0NksXWt4HC8V1vBI9YpMb6h8AAASwtG/ytvapBNOCDvcFTo+/lgaPXogrY2uaw2PRSt3qt3alFnDI5UrMn0hfAAAHBftr/JZE8fE/sUYZdMVN9bu6JKKa3ik86qqdLsAABzX21/lMYkQPPLkl5HV0qXuBI8uhb5czTglPyWCh5Teq6pS+QAAOG7Af5W/8ELHaqU9uFntSAepWJHpDyofAADHDeivcmPCgsd9upngEUWqVWT6g8oHACAh4vqrPMpKpZLU0iLlps/3K3pB+AAAJEy/9zoZMkQ6ejTsMNWOzES3CwAguYwJCx5naIeMrB59lOCRiah8AACS4+23pTPOCDtMtSPzET4AAO5LgbU7kDx0uwAA3BVxUGlARlbHjhE8skHM4WP9+vWaO3euioqKZIzR888/Hzx39OhR3X777SovL9ewYcNUVFSka665Rvv373eyzQCAdHTllb3MZjGyVvJ63W8W3Bdz+GhubtaUKVP04IMPhp1raWnRli1bdOedd2rLli1auXKldu3apb//+793pLEAgDRljPTccyGHvqsnZGRVXU21I9sYa+P/T26M0apVq/T1r3896jU1NTU666yztHfvXpWU9L2hUGNjo3w+n/x+v/Ly8uJtGgAgFTQ1SRH+X87YjswTy/d3wgec+v1+GWM0cuTIRH8UACAG9f5W1TY0q6xgWGJWx2RQKaJIaPj47LPPtHDhQl111VVRU1BbW5va2tqCrxsbGxPZJACAHN7uPpIIwWO4mtSs4QoEouYSZImEzXY5evSorrjiCgUCAT388MNRr6uqqpLP5wv+FBcXJ6pJAABF3+6+3t868F9+//1RB5U2a7isJXggQeHj6NGj+uY3v6na2lqtWbOm176fyspK+f3+4M++ffsS0SQAQCfHtrvvyRjp5ptDDr2or8nI6rXX6GZBN8e7XbqCx3vvvad169YpPz+/1+tzcnKUk5PjdDMAAFEMeLv7nqyVPOF/yzK2A9HEXPk4fPiwtm3bpm3btkmSamtrtW3bNtXV1enYsWO6/PLLtXnzZj3zzDNqb2/XgQMHdODAAR05csTptgMA4jCg7e57MobggZjFPNX2lVde0fnnnx92fP78+brrrrtUVlYW8X3r1q3Teeed1+fvZ6otALij3t8a23b3PUUYvFGu7dqpckJHFkroVNvzzjtPveWVASwbAgBwUb+3u+/p9delc84JO0y1A/3FxnIAgP7rZe2OHTsiblILhCF8AAD6J+q+LFQ7EBt2tQUA9K64OGrwmDyZ4IHYET4AANEZI33wQcih+fqVjKysld56K0ntQlqj2wUAEO5vf5PGjQs7TDcLnEDlAwAQypiowWPfPoIHBo7KBwBkuZDdbUeGr3I6VM1q1VBCBxxD+ACALBayu20goKrJF2ne9jXB80ZWF18srV6dxEYi49DtAgBZKmx3W49Hi2YvUP2IfP1/ujQ4qJTgAadR+QCALFV7sCl8d1uPV2Uj31VbUz7dLEgYKh8AkI2MUdmZp8kTCIQctgGjHa8zvgOJRfgAgGzTuWBYYdPHqlr9gLyBdkkdweMX3zhDEz4Xx34vQAwIHwCQLdasCVupdN72Ndr76EX68tGz9cYd52ve9JIkNQ7ZhDEfAJANIiyPflSDNERHZRslaWDVjpDpuvHslIusQvgAgEyX4A3hQqbrGqmqopwKCnpFtwsAZKpeNoQ7etSZ4BE2XddKi1buVL2/deC/HBmL8AEAmSjChnAVWhFcu2OQQ3Xv2obm8Om61mpPQ4szH4CMRPgAgEzywQdRqx2X/arC8Sm0ZQXD5OnxcV5jVFoQvkw70IXwAQCZwpiOrpaehzurHfPnO/+Rhb5cVVWUy9sZeLzGaHHFGQw6Ra8YcAoAmSBCtWOI2nRUQxK+YNi86SWaNXGM9jS0qLRgKMEDfaLyAQDp7Ec/itrN0hZIfPDoUujL1YxT8gke6BfCBwAkQb2/VRvfbxjYrBBjpH/7t5BD/6Ybgt0sETIJkBLodgEAlw14XYyjR6UhQ8IOG1mtWSPZrzrYWCABCB8A4KJo62LMmjimf10WUcoZXdUOIB3Q7QIALopnXYxgF01eQdi5YtURPJB2qHwAgIu61sU4PoD0ti5GdU2dKldsV0BGnuueUtXqBzRv+xpJVDuQvqh8AICLYlkXo97fqsr/eksBdVwb8Hi0aPYCbR1xWlKChyODZAFR+QAA1/V3XYzayV9S4MqqkGPtHq/+vHK964NK2TwOTqLyAQBJ0Ou6GDk5kjEqO7RfnkAg5JTXGJ0/3d2ly9k8Dk4jfABAKjFGOnJEklTY9LGqVj8gdeaPZC1dzuZxcBrdLgCQCnbvliZMCDt8xfaXtH99a1KXLo91kCzQFyofAJBsxkQMHl2DSpO9dDmbx8FpVD4AIJkiLBrm1THt+9ArW5SE9kTB5nFwEpUPAEiGa66JuiFcu/WqKIWCR5dkV2CQOah8AIDbIoSOf9FPtPOb/yJbnYT2AC4jfACAW1pbpaHhgzRZqRTZhvABAG5gQzggiDEfAJBoEYLHifqbDjcRPJCdCB8AkCjLl0cdVHrQnqjhw5PQJiAF0O0CAIkQIXS8pcn63U/ekv2XJLQHSCGEDwBwWpRqh7XSlCQ0B0g1dLsAwAAFt5rPK+g1eADoQPgAgAGorqnTzCVrddXjmzTzuqdUPfmi4Llz9ZqOHSV4AD0RPgAgTmFbzXs8WjR7gepH5MvIaoM9V4Po3AbCED4AIE61k78UvtW8x6tHq96l2gH0gkwOAPEwRmUj8uUJBBTwdP8d5zVG132LreaB3lD5AIBY3H13cFBpYdPHqlr9gLyBdklsNQ/0F5UPAOivCDNZ/rT9cr326kXa+zFbzQP9RfgAgL60tKh+XIlqS8pVdmi/Cps+Vv2IfJWN+otqtw9ToS9XRSMJHUB/ET4AoDfDhqn68zNVed1TCng88gQC+oe312pl+Vc1Tps0c4lUVVGuedNLkt1SIG0w5gMAojFG9d5cVc6+ITioNODxaEX5V9U1mSVgpUUrd6re35q8dgJphvABAD2tWxcc31E7qihkNksk7dZqT0OLGy0DMgLdLgBwvB6DSssO7e8ob3jCB5t28Rqj0gKm1wL9ReUDACTJ2oizWYqaGrT0G5Pl7TznNUb/+MXPhbxmei0QGyofAHDZZdILL4Qd/mCflR0vSSWaNXGM9jR0T6f959mTQl4D6D/CB4CMUe9vVW1Ds8oKhvU/EESodhSrTvtsscYfd6zQlxvyO3u+BtB/hA8AGaG6pi64yZvH9GP6a22tdPLJYYfv+pnVvrvia0Nc4QfIQoQPAGkvbHfZzumvsyaOiRwCIlQ7qvVNzbPVuivONsQcfoAsxoBTACmh3t+qje83xLVeRm1Dc/justGmv0YIHkYBzbPVMX9ul2jhh7U/gMgIHwCSrrqmTjOXrNVVj2/SzCVrVV1TF9P7ywqGhc2EDZv+umRJxODx6SEra6NPo+2PmMIPAMIHgORyompQ6MtVVUV59OmvxkiVlSHvOVevSdZq5MiB/zv0K/wACGLMB4Ck6q1qEMugzXnTw6fD6vBhacSIsGv//VGrDdcOrN09B5dWVZRr0cqdareWtT+APhA+ACRVV9Xg+AASb9UgZPprhC6Wj1SgMfYjDTB3RB1cGhZ+AEQUc7fL+vXrNXfuXBUVFckYo+effz7kvLVWd911l4qKipSbm6vzzjtPb7/9tlPtBZBh+uwyiUeE4JGjzzTGfhT/7+zUWzdRoS9XM07JJ3gAfYg5fDQ3N2vKlCl68MEHI57/xS9+ofvuu08PPvigampqNG7cOF100UVqamoacGMBZKZ500u0YeH5+u33ztaGhefHP0X1D3+IGDw+a7VqszkDbGUHBpcCAxdzt8ucOXM0Z86ciOestbr//vt1xx13qKKiQpK0fPlyjR07Vs8++6yuvXagxU4AmWrAK4ZGCB3f17/rMft9nTCAdvXkZDcRkK0cne1SW1urAwcO6OKLLw4ey8nJ0Ve+8hVt3Lgx4nva2trU2NgY8gMA/RYIRAwe/+/zVo/Z7zv+cQnpJgKyjKMDTg8cOCBJGjt2bMjxsWPHau/evRHfU1VVpZ///OdONgNAljhWXKpBH0T4f4u1uiyBn8vgUmBgErLOh+nxV4i1NuxYl8rKSvn9/uDPvn37EtEkAJnGmLDgMV77JGujvMFZDC4F4udo5WPcuHGSOioghYWFweMHDx4Mq4Z0ycnJUU6OMwPBAGSB7dulKVPCDrcfs/rAm4T2AIiZo5WPsrIyjRs3TmvWrAkeO3LkiF599VWdc845Tn4UgGxkTFjw+L0ukayVl+ABpI2YKx+HDx/W7t27g69ra2u1bds2jR49WiUlJbrpppu0ePFiTZgwQRMmTNDixYs1dOhQXXXVVY42HECWidB1+/Iaq0u/moS2ABiQmMPH5s2bdf755wdf33LLLZKk+fPn61e/+pV+/OMfq7W1VT/84Q916NAhfelLX9JLL72kERGWOAaAvhy6/HsateKJ8BPWitwBpCdjrUujs/qpsbFRPp9Pfr9feXl5yW4OgGSKUO34ZtFr+s8Pz01CYwD0Jpbvb/Z2AZB6PvlEys8PO2wDVv8ZeeIcgDRC+ACQWqJMy5e1IncAmSEh63wAQFwiBI+tm464tnYHAHcQPgAk3a6F/xG54mGtpp412P0GAUgowgeA5DJGE5d+N+TQQ6c9RLUDyGCM+QCQFIGj7fIMifC/IGt1vfvNAeAiwgcA9xkTuexKtQPICnS7AHBXhLEddTV/I3gAWYTwASBm9f5WbXy/QfX+1n6/58XK9VEHlZZMO9HB1gFIdXS7AIhJdU2dKlfuUMBKHiNVVZRr3vSS3t9kjOb0OLT9xAs1+W8vJ6ydAFIXlQ8A/Vbvbw0GD0kKWGnRyp1RKyCffaao1Q6CB5C9CB8A+q22oTkYPLq0W6s9DS1h1643s3RCbuTgASC7ET4A9FtZwTB5euQJrzEqLRgaetAYzdJrIYcO/d8NBA8AkggfAGJQ6MtVVUW5vJ1dKV5jtLjiDBX6ciVJ99/2YdRullH/a6abTQWQwoy1qfWnSCxb8gJIjnp/q/Y0tKi0YGgwePS2IRyAzBfL9zezXQDErNCXGwwdhw5Jo0ZHCB7Hjkler8stA5AO6HYBELefmrsjBw9r+wwe8awVAiAzUPkAEB9jdHePQ233PqCcWxf0+da41goBkDGofACISeWtR6IOKu1P8Ih1rRAAmYfKB4D+M0ZVkY7HMKi0t7VCgoNXAWQ0Kh8A+vThh4pc7WhoiHk2S7/XCgGQsQgfAHp1nnlFnxsfZVBpfn7Mv6+vtUIAZD66XQBEZ4xe6XEocN0P5Hnk4QH92nnTSzRr4pjwtUIAZAXCB4Aw3/iG9F+/i1ztcKpcevxaIQCyC+EDQIjVZrb+Sy+Fn2ClUgAOYcwHAEnSn/8syRjN7hk8du4keABwFJUPACoxdarTSeEnCB0AEoDKB5DtjAkLHnbCBIIHgIQhfABZaupURV67IxCQ2bXL9fYAyB6EDyAL3WHu0dZtUdbuiBRIAMBBhA8gi7z2miRjdI9+Enrid7+jmwWAaxhwCmSJXNOqVkVYwpzQAcBlhA8gwwUCksdrFHHPWIIHgCSg2wXIYMOGdQSPME1NBA8ASUP4ADLUHPOimluiDCodPtz9BgFAJ8IHkGFWrJBkjF7UJaEnbr+dageAlMCYDyCDGGMVces3QgeAFELlA8gAbW3S/5izCB4A0gLhA0hzxkg5JxidpZrQE7t2ETwApCS6XYA09nmzW1YTwk8QOgCkMCofQBp66CFJxmh3z+Bx1lkEDwApj8oHkGaMkawibwjHviwA0gGVDyCJ6v2t2vh+g+r9EdcfDeH3S0vMwsjBgw3hAKQRKh9AklTX1Kly5Q4FrOQxUlVFueZNL5HUEUpqG5pVVjBMhb7cYLVjYc9f8uKL0te+5nrbAWAgCB9AEtT7W4PBQ5ICVlq0cqdmTRyj9bs+Cgkl/j9MlNXE8F/C2A4AaYpuFyAJahuag8GjS7u12rL3UFgoGTn7z6ofkR96McEDQBojfABJUFYwTJ4eQzS8xihgbXgo8Xi1Z2SRJKn+wCfauPujfo0RAYBURfgAkqDQl6uqinJ5OweJeo3R4oozVHzCaPVMH95Au0o/3a/q/9mrmfdv1FWPb9LMJWtVXVOXjKYDwIAx5gNIknnTSzRr4hjtaWhRacFQFY3MlZXR0skXadHsBWr3eOUNtGvx0P3Svn2qXLI24hiRQl9ucv9FACBGhA8giQp9uZ2zWWxwCu287Ws0q3aL9owsUumOTSr05Wrj+w0Rx4jsaWghfABIO3S7AEk0b570npkQtiFcYdPHmlG3PRgsoo0RKS0Y6lZTAcAxhA8gSYyRqv/TaIJ2h57YuzdsNku0MSJUPQCkI7pdAJf95S/S5V/YKavy8JO9TKHtOUaE4AEgXRE+ABd1rVS6s+eJ2bOlP/6xz/d3jREBgHRG+ABcEnVDOBYMA5BlGPMBJNi0adLD5ocEDwDoROUDSKCo1Y5XXpG+8hXX2wMAqYDKB5AAGzZIo8yh6NUOggeALEblA3BYV7XjUKSTdLMAAJUPwCmBQC/dLG1tBA8A6ET4ABwwfLj0v71PR+9mGTLE/UYBQIoifAADZIx0uNnoac0PPbFsGdUOAIjA8fBx7Ngx/eQnP1FZWZlyc3N18skn6+6771YgEHD6o4CkWrVK8pr26NWOG290v1EAkAYcH3C6dOlSPfroo1q+fLlOP/10bd68Wd/+9rfl8/n0ox/9yOmPA5LCGKlB+WrXJ+EnqXYAQK8cDx+vv/66LrvsMl166aWSpNLSUv32t7/V5s2bnf4owHVHjkg5OVEGlR44II0d636jACDNON7tcu655+pPf/qTdu3aJUl66623tGHDBl1yySURr29ra1NjY2PID5CKjJHOydkcvZuF4AEA/eJ45eP222+X3+/XqaeeKq/Xq/b2dt1zzz268sorI15fVVWln//85043A3BU1Cm03/ymVF3tfoMAII05Xvmorq7Wb37zGz377LPasmWLli9frnvvvVfLly+PeH1lZaX8fn/wZ9++fU43CYjbk0/2sSEcwQMAYuZ45eO2227TwoULdcUVV0iSysvLtXfvXlVVVWn+/Plh1+fk5CgnJ8fpZgADZoz0qK6V1WPhJxlUCgBxczx8tLS0yOMJLah4vV6m2iJtNDVJeXlRqh2bNklnneV+owAggzgePubOnat77rlHJSUlOv3007V161bdd999+s53vuP0RwGOM0bKV4OsxoSfpNoBAI4w1jr7f9SmpibdeeedWrVqlQ4ePKiioiJdeeWV+ulPf6oh/VhiurGxUT6fT36/X3l5eU42DehV1LEdJSXS3r3uNwgA0kgs39+Oh4+BInzAbXffLf3sZ1GCx7FjktfrfqMAIM3E8v3N3i7IasZIe3/2ZPTZLAQPAHAc4QNZ6cCB7m6WJ/V/Qk/+5jeM7wCABHJ8wCmQ6oyRvDomq8HhJwkdAJBwVD6QVbo2hDtG8ACApCF8ICt873vd3Sz5PXei/eSTfgWPen+rNr7foHp/a4JaCQDZgW4XZDxjpC/pDVnNCD/Zz2pHdU2dKlfuUMBKHiNVVZRr3vQSh1sKANmBygcy1nvvdVc73ugZPG66qd/Bo97fGgwekhSw0qKVO6mAAECcqHwgI5nOmbNRp9DGoLahORg8urRbqz0NLSr05cbZQgDIXlQ+kHGMkZ7VlY4ED0kqKxgmT49f5TVGpQVD42whAGQ3wgcyxgUXdHezXKnnQk/u3Bn3bJZCX66qKsrl7SyneI3R4oozqHoAQJzodkFGMEYq1H5ZfS78pANTaOdNL9GsiWO0p6FFpQVDCR4AMABUPpDWNm3qrnbs7xk8pkxxdO2OQl+uZpyST/AAgAGi8oG01eug0kCg+wIAQEqh8oG0Y21HrrhB/xZ9UCnBAwBSFpUPpJXSUmnv3ijVjlWrpK9/3e0mAQBiRPhA2jBGGqI2WZ0QfpJ9WQAgbdDtgpT3+993DyptI3gAQNqj8oGU1uug0qYmafhwdxsEABgwKh9ISUeOdASP87Qu+qBSggcApCXCB1KOMVJOTke1Y50uCD155510swBAmqPbBSmlo5vFykbKxYQOAMgIVD6QEp56qiN4/F6XEDwAIMNR+UDS9TqodPdu6ZRT3G0QACChqHwgaZqaOoLHSdoTfVApwQMAMg7hA66p97dq4/sNqve3yhgpL6+j2rFHZaEXnnce3SwAkMHodoErqmvqVLlyhwJWsgFp+ORyNW0/KfxCNoQDgIxH5QMJV+9vDQYPSTIe6cTZ21Q/Ij/0QjaEA4CsQPhAwtU2NAeDR5d2j1d7RhZ1vFi9mm4WAMgidLsgoT76SJo1dZg+d11HxaOLN9Cu0k/3EzoAIAtR+UDCGCOdeKLU2pSnX6xeJm+gXVJH8Fi8+kEVNjYkuYUAgGQw1qbWn56NjY3y+Xzy+/3Ky8tLdnMQp0hrd9SPyNeekUUqrXlVhWNHJallAIBEiOX7m8oHHHXzzR3B40K9HLZ2R2HTx5pRt53gAQBZjjEfcEyvK5U+8IC0YIG7DQIApCTCBwbsr3/tWoiUDeEAAH2j2wUDYkxH8GBDOABAf1H5QNx67Wb54APpc59zt0EAgLRA5QMxu/zyjuBxst6PviEcwQMAEAWVD8Sk12rH1VdLTz/tboMAAGmH8IF+2bJFOvPMjn+OWu0AAKAf6HZBn4zpCB4/0b8QPAAAA0blA1FZK3k642nE0PH669LZZ7vbKABA2qPygYi++MWO4DFcTdGrHQQPAEAcqHwgTK+DSsePl/btc7dBAICMQuUDQS+/3EfwOHKE4AEAGDDCByR1hI6LLpLm6oXo3SyDB7vfMABAxiF8ZKF6f6s2vt+gen+rjh0LrXa8oMtCL/6P/2A2CwDAUYz5yDLVNXWqXLlDASvJSh//sVxG4xWQN/xiQgcAIAGofGSRen9rd/CQJCMVzH5LH444MfxiggcAIEEIH1mktqG5O3h0sh6P9ows6j5w8CDBAwCQUISPDHf8+I5ZU4fJBkLPewPtKv10f8cLa6UxY9xvJAAgqxA+Mlh1TZ1mLlmrqx7fpLPvWavcso/0i9XL5A20S+oIHotXP6jCH99EtQMA4BpjbWp96zQ2Nsrn88nv9ysvLy/ZzUlb9f5WzVyyNqSbxRto14ZHvyNJ2jOySKWf7ldhY0Pw+tqGZpUVDFOhLzcZTQYApLFYvr+Z7ZKhIo3vaPd4tWdkkWbs26HCpo+D1Y7jZ8B4jFRVUa5500uS0GoAQDag2yVBjh9r4bbly9X7+I4//zkYPHrOgAlYadHKnUlpNwAgO1D5SIBkVhK6FgzL0xH9YvUyLZq9QO0eb/f4js5uli4RKyTWak9DC90vAICEIHw4LFolYdbEMQn9Mm9slHy+jn8+pJEaKb+0XZpVu6VjfMdZ5Sp866Ww95UVDJPHKHRsiDEqLRiasLYCALIb3S4O662SkCiTJnUHDyvTETw6FTZ9rBm1W1X4u2civrfQl6uqinJ5O0smXmO0uOIMqh4AgISh8uEwtysJXd0sM7RRGzUz/IJ+TGaaN71EsyaO0Z6GFpUWDCV4AAASisqHw9yqJCxfHrohXFjwePHFmNbuKPTlasYp+QQPAEDCUflIgERXErpCh0ftao/0nzC1lm4BACAElY8ESUQl4W9/6w4ev9Qt4cHjggsIHgCAlEflI03k5kqffdbxz1Ym/ILDh6Vhw9xtFAAAcUhI5ePDDz/Ut771LeXn52vo0KH6u7/7O7355puJ+KisYExH8ChVbeTgYS3BAwCQNhwPH4cOHdLMmTM1ePBgvfjii3rnnXf0y1/+UiNHjnT6ozLeL38ZOqi0VieHXvDMM2HdLMlcWRUAgP5wvNtl6dKlKi4u1lNPPRU8Vlpa6vTHZDxzXIEjarWjB/ZoAQCkA8crHy+88IKmTZumb3zjGzrxxBM1depUPf7441Gvb2trU2NjY8hPNqut7Q4e39UT4cEjNzdi8GCPFgBAunA8fPz1r3/VI488ogkTJmj16tW67rrrdOONN+rpp5+OeH1VVZV8Pl/wp7i42OkmpQ1jpJM7e1asjJ7Q90Iv+OADqSXySqnJWFkVAIB4GGudnZs5ZMgQTZs2TRs3bgweu/HGG1VTU6PXX3897Pq2tja1tbUFXzc2Nqq4uFh+v195eXlONi1lWSt5OmNgnvzya2Tki3pR72/VzCVrw1ZW3bDw/LDpvvX+VtU2NKusYBiLigEAHNHY2Cifz9ev72/HKx+FhYU67bTTQo594QtfUF1dXcTrc3JylJeXF/KTTW6/vTt41GhaePC4885+rd3R35VVq2vqNHPJWl31+CbNXLJW1TWR/7sAAJAojg84nTlzpt59992QY7t27dJJJ53k9EelvT4HlQYCoRf1oa+VVZO14y4AAMdzvPJx880364033tDixYu1e/duPfvss3rsscd0/fXXO/1RaWv79u5McZ7WRZ/NEkPw6NLbyqqMCwEApALHKx/Tp0/XqlWrVFlZqbvvvltlZWW6//779U//9E9Of1Ra6rPasXmzdOaZCflst3fcBQAgEscHnA5ULANW0kl7uzSoM+p5dUzHNDj8Ihf+U1TX1GnRyp1qtzY4LoS1QAAAAxXL9zd7u7jgmmukX/+6458f1bW6Vo+FXjB3rvTCC660JdE77gIA0BfCR4L12c3S0tKxcJiLCn25hA4AQNIkZGM5SBs2dAePCdoVfVCpy8EDAIBko/KRAH1WO1askCoq3GsQAAAphPDhoCNHpJyc7tf93RAOAIBsQreLQ269tTt4/INWhgePsWMJHgAAiMqHI/rsZjlwoCN8AAAAKh8DsXFjd/AYocbo3SwEDwAAgggfx6n3t2rj+w2q97f2ea0x0syZHf+8XNeoUb7QC555hm4WAAAioNulU3VNXXDTNY+RqirKI6782dwsDR/e/dqJDeEAAMgmVD4UfbfXnhWQq6/uDh5naVN48Dj55Lg3hAMAIFtQ+VDvu712rQTa56DS99/vCB8AAKBXVD7Uvdvr8bp2e129ujt4eHUs+qBSggcAAP1C+FDHXidVFeXydqaMrt1ei0bm6mtf67hmke4J34n27rsZVAoAQIzodul0/G6vowYN1RdKu/dciVjtOHJEGjw4/DgAAOgVlY/jFPpyddcP84PBo1S10btZCB4AAMSFysdxjh9UuksTNEG7Qy94/XXp7LPdbRQAABmG8KGOlUq7FgyTrGykghBjOwAAcETWd7vk53cHjyv02/DgcdVVBA8AAByUtZWPQ4ek0aO7X3+mHOXoSOhFfr+Ul+duwwAAyHBZWfn453/uDh5D1SwrEx48rCV4AACQAFkVPrpWPv/lLzte36p71azhoRe9/DLdLAAAJFBWdbt4jotaUafQAgCAhMqaykdXrjhVfw4PHrfeSvAAAMAlWVP5MEbyL1qqvMULQ0989JFUUJCcRgEAkIWyJnzI2vDgQbUDAADXZU23i4yRbrml45//+EeCBwAASWKsTa1v4cbGRvl8Pvn9fuUx1RUAgLQQy/d39lQ+AABASsieMR+S6v2t2rznExljdOZJo1Toy012kwAAyDpZEz6qa+q0cMUOdfUxGUlL/rFc86aXJLNZAABknazodqn3t4YED0mykipX7lC9vzVZzQIAICtlRfiobWhWpFG1ASvtaWhxvT0AAGSzrAgfZQXDIi2mLo+RSguGut4eAACyWVaEj0Jfrpb8Y3lIADFGqqooZ9ApAAAuy5oBp/Oml2jWxDF6c88hGSN9kdkuAAAkRdaED6mjAvK/phA4AABIpqzodgEAAKmD8AEAAFxF+AAAAK4ifAAAAFcRPgAAgKsIHwAAwFWEDwAA4CrCBwAAcBXhAwAAuIrwAQAAXEX4AAAArkq5vV2stZKkxsbGJLcEAAD0V9f3dtf3eG9SLnw0NTVJkoqLi5PcEgAAEKumpib5fL5erzG2PxHFRYFAQPv379eIESNkjIl4TWNjo4qLi7Vv3z7l5eW53ML0xX2LH/cuPty3+HDf4se9i48T981aq6amJhUVFcnj6X1UR8pVPjwej8aPH9+va/Py8ni44sB9ix/3Lj7ct/hw3+LHvYvPQO9bXxWPLgw4BQAAriJ8AAAAV6Vl+MjJydHPfvYz5eTkJLspaYX7Fj/uXXy4b/HhvsWPexcft+9byg04BQAAmS0tKx8AACB9ET4AAICrCB8AAMBVhA8AAOCqlA0fDz/8sMrKynTCCSfozDPP1Guvvdbr9a+++qrOPPNMnXDCCTr55JP16KOPutTS1BLLfXvllVdkjAn7+ctf/uJii5Nv/fr1mjt3roqKimSM0fPPP9/ne3jeOsR673jmpKqqKk2fPl0jRozQiSeeqK9//et69913+3wfz1x8945nTnrkkUc0efLk4AJiM2bM0IsvvtjrexL9vKVk+KiurtZNN92kO+64Q1u3btWXv/xlzZkzR3V1dRGvr62t1SWXXKIvf/nL2rp1qxYtWqQbb7xRK1ascLnlyRXrfevy7rvvqr6+PvgzYcIEl1qcGpqbmzVlyhQ9+OCD/bqe561brPeuSzY/c6+++qquv/56vfHGG1qzZo2OHTumiy++WM3NzVHfwzPXIZ571yWbn7nx48dryZIl2rx5szZv3qwLLrhAl112md5+++2I17vyvNkUdNZZZ9nrrrsu5Nipp55qFy5cGPH6H//4x/bUU08NOXbttdfas88+O2FtTEWx3rd169ZZSfbQoUMutC49SLKrVq3q9Rqet8j6c+945sIdPHjQSrKvvvpq1Gt45iLrz73jmYts1KhR9oknnoh4zo3nLeUqH0eOHNGbb76piy++OOT4xRdfrI0bN0Z8z+uvvx52/ezZs7V582YdPXo0YW1NJfHcty5Tp05VYWGhLrzwQq1bty6RzcwIPG8DxzPXze/3S5JGjx4d9Rqeucj6c++68Mx1aG9v13PPPafm5mbNmDEj4jVuPG8pFz4aGhrU3t6usWPHhhwfO3asDhw4EPE9Bw4ciHj9sWPH1NDQkLC2ppJ47lthYaEee+wxrVixQitXrtSkSZN04YUXav369W40OW3xvMWPZy6UtVa33HKLzj33XJ1xxhlRr+OZC9ffe8cz12HHjh0aPny4cnJydN1112nVqlU67bTTIl7rxvOWcrvadjHGhLy21oYd6+v6SMczXSz3bdKkSZo0aVLw9YwZM7Rv3z7de++9mjVrVkLbme543uLDMxdqwYIF2r59uzZs2NDntTxzofp773jmOkyaNEnbtm3Tp59+qhUrVmj+/Pl69dVXowaQRD9vKVf5KCgokNfrDftr/eDBg2FJrMu4ceMiXj9o0CDl5+cnrK2pJJ77FsnZZ5+t9957z+nmZRSeN2dl6zN3ww036IUXXtC6des0fvz4Xq/lmQsVy72LJBufuSFDhujzn/+8pk2bpqqqKk2ZMkXLli2LeK0bz1vKhY8hQ4bozDPP1Jo1a0KOr1mzRuecc07E98yYMSPs+pdeeknTpk3T4MGDE9bWVBLPfYtk69atKiwsdLp5GYXnzVnZ9sxZa7VgwQKtXLlSa9euVVlZWZ/v4ZnrEM+9iyTbnrlIrLVqa2uLeM6V582xoasOeu655+zgwYPtk08+ad955x1700032WHDhtk9e/ZYa61duHChvfrqq4PX//Wvf7VDhw61N998s33nnXfsk08+aQcPHmx/97vfJetfISlivW//+q//aletWmV37dpld+7caRcuXGgl2RUrViTrXyEpmpqa7NatW+3WrVutJHvffffZrVu32r1791pred56E+u945mz9gc/+IH1+Xz2lVdesfX19cGflpaW4DU8c5HFc+945qytrKy069evt7W1tXb79u120aJF1uPx2Jdeeslam5znLSXDh7XWPvTQQ/akk06yQ4YMsV/84hdDplLNnz/ffuUrXwm5/pVXXrFTp061Q4YMsaWlpfaRRx5xucWpIZb7tnTpUnvKKafYE044wY4aNcqee+659ve//30SWp1cXVPxev7Mnz/fWsvz1ptY7x3PnI14vyTZp556KngNz1xk8dw7njlrv/Od7wS/F8aMGWMvvPDCYPCwNjnPm7G2cxQJAACAC1JuzAcAAMhshA8AAOAqwgcAAHAV4QMAALiK8AEAAFxF+AAAAK4ifAAAAFcRPgAAgKsIHwAAwFWEDwAA4CrCBwAAcBXhAwAAuOr/B4MYmxhx04/ZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Your code here\n",
    "#A=np.array([[1,2,3],[4,5,6],[7,8,9]])\n",
    "#x = A[0, :] #从一个矩阵中提取出一行作为一个向量\n",
    "#y1 = np.array([2, 3, 5])\n",
    "#plt.plot(x, y1) #画出折线图\n",
    "#y2 = np.array([2.5, 2.8, 5.3])\n",
    "#plt.plot(x, y2, '.') #画出散点图\n",
    "#plt.show()\n",
    "\n",
    "test_x = test[:, 0]\n",
    "test_y = test[:, 1]\n",
    "\n",
    "train_y_1 = omega_1 * test_x + b_1\n",
    "train_y_2 = omega_2 * test_x + b_2\n",
    "\n",
    "plt.plot(test_x, train_y_1, 'b')\n",
    "plt.plot(test_x, train_y_2, 'r')\n",
    "plt.plot(test_x, test_y, '.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f227579c",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">4) 在训练数据集'train2.csv'上求一个三元线性回归模型$\\hat{y}=w_0 + w_1 x_1 + w_2 x_2 + w_3 x_3$的使得损失函数$l(w_0,w_1,w_2,w_3)=\\frac{1}{2}\\sum_{i=1}^m(\\hat{y}^{(i)}-y^{(i)})^2$最小的参数$w_0,w_1,w_2$以及$w_3$。并在测试数据集'test2.csv'上进行预测，输出预测结果的均方误差$MSE(\\hat{y},y)=\\frac{1}{n}\\sum^n_{i=1}(y^{(i)}-\\hat{y}^{(i)})^2$, $n$为测试集中样本个数。</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a1f6178",
   "metadata": {},
   "source": [
    "方法① 同2)中的方法③。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4ea8a93",
   "metadata": {},
   "source": [
    "方法② 类似2)中的方法②。算法步骤如下：1.初始化模型参数$w_0,w_1,w_2,w_3$的值；2.在负梯度的方向上更新参数(批量梯度下降、小批量随机梯度下降或者随机梯度下降均可)，并不断迭代这一步骤，更新公式(以小批量随机梯度下降为例)可以写成：$$w_j\\gets w_j-\\frac{\\eta}{\\left|B\\right|}\\sum_{i\\in{B}}x_j^{(i)}(w_0 + w_1 x_1^{(i)}+w_2 x_2^{(i)}+w_3 x_3^{(i)}-y^{(i)}), j=0,1,2,3$$, 其中$x_0^{(i)}=1$， 其中$\\eta$表示学习率,$B$表示每次迭代中随机抽样的小批量，$\\left|B\\right|$则表示$B$中的样本数量。3. 终止条件为迭代次数达到某一上限或者参数更新的幅度小于某个阈值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "115f58e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "omega = [4.483974672173968, 1.2308467502840124, 2.17258322878662, 3.20945679025689]\n",
      "MSE = 0.383779830817784\n"
     ]
    }
   ],
   "source": [
    "# Your code here\n",
    "\n",
    "train_frame = pd.read_csv('train2.csv')\n",
    "test_frame = pd.read_csv('test2.csv')\n",
    "\n",
    "train = np.array(train_frame)\n",
    "test = np.array(test_frame)\n",
    "\n",
    "# 将数据集分成k个mini_batch\n",
    "def get_mini_batches(train, k):\n",
    "    np.random.shuffle(train)\n",
    "    mini_batches = [train[i:i + k] for i in range(0, len(train), k)]\n",
    "    return mini_batches\n",
    "\n",
    "k = 10\n",
    "mini_batches = get_mini_batches(train, k)\n",
    "\n",
    "# 学习率\n",
    "learning_rate = 0.01\n",
    "# 阈值 >= 0.001\n",
    "threshold = 0.001\n",
    "\n",
    "def gradient_descent(omega, mini_batches):\n",
    "    # omega的更新公式\n",
    "    def omega_gd(omega, mini_batch):\n",
    "        res = omega.copy()\n",
    "        mini_batch_with_x0 = np.insert(mini_batch, 0, 1, axis=1)\n",
    "        total = 0\n",
    "        for i in range(len(res)):\n",
    "            for row in mini_batch_with_x0:\n",
    "                xj = row[i]\n",
    "                total += xj * (omega[0] + omega[1] * row[1] + omega[2] * row[2] + omega[3] * row[3] - row[4])\n",
    "            res[i] -= learning_rate * total / len(mini_batch)\n",
    "        return res\n",
    "\n",
    "    while True:\n",
    "        random_index = np.random.randint(0, len(mini_batches))\n",
    "        mini_batch = mini_batches[random_index]\n",
    "        omega_new = omega_gd(omega, mini_batch)\n",
    "        flag = True\n",
    "        # 终止条件为参数更新的幅度小于阈值threshold\n",
    "        for x, y in zip(omega_new, omega):\n",
    "            if abs(x - y) >= threshold:\n",
    "                flag = False\n",
    "                break\n",
    "        omega = omega_new\n",
    "        if flag:\n",
    "            break\n",
    "    return omega\n",
    "\n",
    "#omega任意初始值\n",
    "omega = [1, 1, 1, 1]\n",
    "omega = gradient_descent(omega, mini_batches)\n",
    "\n",
    "test_x = test[:, :3]\n",
    "train_y = [omega[0] + omega[1] * x1 + omega[2] * x2 + omega[3] * x3 for x1, x2, x3 in test_x]\n",
    "test_y = test[:, 3]\n",
    "MSE = np.sum((test_y - train_y) ** 2) / len(test_y)\n",
    "print('omega = ' + str(omega))\n",
    "print('MSE = ' + str(MSE))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1507673b",
   "metadata": {},
   "source": [
    "<font color=blue size=4>2.相关链接</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6892053c",
   "metadata": {},
   "source": [
    "1.实验报告提交链接(有效期直至9.15 14:20): https://send2me.cn/211f55kq/Sgav4JPN0foh9Q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ad2992c",
   "metadata": {},
   "source": [
    "2.实验课件获取链接: https://www.jianguoyun.com/p/DWcqLm4Qp5WhChjZi5sFIAA"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
